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Abstract. Calculations of the ground state of inhomogeneous many-electron systems involve 
a solving of the Poisson equation for Coulomb potential and the Schrodinger equation for single- 
particle orbitals. Due to nonlinearity and complexity this set of equations, one believes in the 
iterative method for the solution that should consist in consecutive improvement of the potential 
and the electron density until the self-consistency is attained. Though this approach exists for a 
long time there are two grave problems accompanying its implementation to infinitely extended 
systems. 

The first of them is related with the Poisson equation and lies in possible incompatibility 
of the boundary conditions for the potential with the electron density distribution renewed by 
means of the Schrodinger equation. Rigorously speaking, it means the iteration process cannot 
be continued. The resolution of this difficulty is proposed for both infinite conducting systems 
in jellium approximation and periodic structures. It provides the existence of self-consistent 
solution for the potential at every iteration step due to realization of a screening effect. 

The second problem results from the existence of continuous spectrum of Hamiltonian 
eigenvalues for unbounded systems. It needs to have a definition of Hilbert space basis with 
eigenfunctions of continuous spectrum as elements, which would be convenient in numerical 
applications. It is proposed to insert a limiting transition into the definition of scalar product 
specifying the Hilbert space. It provides self-adjointness of Hamiltonian and, respectively, the 
orthogonality of eigenfunctions corresponding to the different eigenvalues. In addition, it allows 
to normalize them effectively to delta- function and to prove the orthogonality of the 'right' and 
'left' eigenfunctions belonging to twofold degenerate eigenvalues. 



1. Problem definition 

Ground-state calculations of inhomogeneous many-electron systems involve generally a solving of 
the Poisson equation for averaged Coulomb potential u(r) at given spatial electron density n(r) 
and the Schrodinger equation for single-particle orbitals V'E(r) in the potential u e ff accounting 
by some approximation for the difference between u(r) and the microscopic local field. In the 
density functional theory the corresponding set of Kohn-Sham equations for the spin-unpolarized 
electron gas has the form (in atomic units \e\ = m = h = 1) 
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Here E is the energy eigenvalue of the single-particle Hamiltonian, Ep is the Fermi energy of 
the electrons, N+ is the density of the positive background, u is the Coulomb potential energy 
of the electron, and u xc is the exchange-correlation potential energy, which is assumed in the 
local-density approximation u xc (r) = u xc [n(r)]. 

Due to nonlinearity and complexity of this set of equations, one believes in the iterative 
solution procedure that should consist in consecutive improvement of u(r) and n(r) until the 
self-consistency is attained. While this approach exists for a long time and as if were used in 
many articles there are two grave problems accompanying its applications and pertinent to very 
background of the iterative method in the case of infinitely extended many-electron systems. 

The first of them is related with Poisson equation and lies in possible incompatibility of the 
boundary conditions for it(r) with the distribution n(r) renewed by means of the Schrodinger 
equation solutions and substituted to the right-hand side of the Poisson equation (J3J) . Rigorously 
speaking, it means the iteration process cannot be continued. To manage this difficulty there 
were several empirical techniques suggested (see e.c. jTj) but their shortcomings are either lack of 
iteration convergence and a transfer to some kind of variation solution instead (see j2j, Appendix 
B), or an appearance of solution instability 3 or a change in the positive charge distribution 
that violates the derivation conditions of self-consistent field equations as variational equations 
. In the case of extended but finite systems, the effect may result in the non-physical growth 
of the electric filed far away from the inhomogeneity region [3]. The compliance of obtained 
solution with true one is an open question in all cases. The suggested approach to deal with 
this problem are described in Section [2] for systems with Fermi level given at the infinity and in 
Section for systems with given number of electrons. 

The second problem results from the existence of continuous spectrum of Hamiltonian 
eigenvalues for unbounded systems. It is necessary to have such a definition of Hilbert space 
with eigenfunctions of continuous spectrum as elements, which would be convenient in numerical 
applications. A limiting transition to Hamiltonian operator with continuous spectrum by means 
of the unlimited increase in system size is used to build up mathematically the Hilbert space of 
singular self-adjoined operators (see e.c. [5]). However, this way is practically inappropriate for 
inhomogeneous systems. 

It is suggested here to introduce a limiting transition into the definition of the scalar product 
specifying the Hilbert space. The adopted form of the scalar product provides self-adjointness 
of Hamiltonian operator and, respectively, the orthogonality of all eigenfunctions corresponding 
to the different eigenvalues. In addition, it allows to normalize them effectively to delta-function 
and to prove the orthogonality of the "right" and "left" current-carried eigenfunctions belonging 
to twofold degenerate eigenvalues. This is particularly of the essence for the problem of tunneling 
through a self-consistent barrier considered as an example in Section |31 Also this approach is 
implemented to the Bloch wave functions of periodic solids in Section |1J 

2. Poisson equation and iteration algorithm 

Let us illustrate the neutrality problem with one-dimensional Poisson equation given on semi- 
axis z £ [0, oo): 
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The simple integration results in 
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The finiteness condition of u(oo) < oo requires \im z ^ 00 zu'(z) = and 
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47T / dz\Z\p(z\) = u(0) — u(oo) (9) 
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Assuming an energy scale such that u(oo) = it is easily to see the existence of the direct 
relationship between the total charge and a boundary condition for the electric field u'(0) in 
Eq.Q or between the total dipole moment and a boundary condition for the potential u(0) 
in Eq. ([?])■ However, the electron density n^(r) in Eq.Q obtained after the solution of the 
Schrodinger Eq.Q at an i-th iteration step may not satisfy the imposed boundary conditions 
for the Poisson equation as usually it takes place. In this case there is no possibility to solve the 
Eq.JSJ) and the iteration procedure should be stopped. 

In the presented approach this difficulty is removed by partition of full density n(r) in two 
terms 

n(z) = n ind [u{z)\ + n Q (z), (10) 

where n; n( j is defined as a function of the Coulomb potential u via its relation to the effective 
potential ii e g- = u(z) + ii xc [n(z)] by the known quasi-classical expression 

2 3/2 3/2 
n-md (u, n) = [E F - u cfi (z, n(z))) 1 (11) 

and 

n Q (z) = n(z) - n ind (z) (12) 

is named as the quantum correction. Using the definition ()1U|) . the Poisson equation © can be 
rewritten in the form 

u" + 47rn ind [u,n Q ] = 4vr(iV + - n Q (z)). (13) 

If the pair functions n{z) and u(z) is the true self-consistent solution of the problem (JIJ1J) then 
the Ea.H13|) is the simply rearranged Eq.©. However, the Ea. (|13l) is much more appropriate 

(i) 

for the iterative procedure since, due to the screening effect, the induced electron density nj nd 
depending on the unknown Coulomb potential provides the existence of self-consistent solution 
«W(r) at every iteration step and any possible spatial dependence of right-hand side of Eq. (|13[) . 
It is useful to note that the expression (fTT|) is a good approximation of the solution of Schrodinger 
equation for the smooth part of the u e s and produces the true screening of the long-range part 
of the Coulomb potential since the riind-is found simultaneously with u in the course of self- 
consistent solution of the Poisson equation. The rest short-range variations of the density n{z) 
are exactly described by nq(z), which is to be found in the usual iterative cycle after the solution 
of Schrodinger equation. This is the root cause of the algorithm efficiency. It will be convenient 
to designate the Eq. (fT3)) as self-consistent Poisson equation. 

The full iteration algorithm in the case of the one-dimensional inhomogeneity of the charge 
distribution can be described now as following 

i = 0,l,...; ng } =0, (14) 
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Here rtg (z) is the electron density that is self-consistent with the Coulomb potential u^\z) at 



the given quantum correction density nq (z)). The self-consistent Poisson equation <|15|) is solved 
as the boundary problem and the Schrodinger equation l|17|l is solved as the Cauchy problem. 

It needs to say a few words relative to the speciality of the quasi-classical expression (jlljl 
for n in( j in the case of density functional approach. Because -u xc depends itself on the electron 
density the Eq. ()ll|) determines n ia & as the implicit function of the Coulomb potential u when 
riQ is known. This function has the physical meaning only under condition dn m &/du < 0, which 
is the stability condition for solutions of the Ea. H15|) (see a discussion in jSj). 

The validity of the described method for semi-infinite electron systems was examined by 
calculations of surface properties of simple metals and quantum corrections to the capacity of 
barrier structures [J). The expected convergency of the iteration procedure has been obtained 
and the true self-consistency of the solution has been successfully checked by the Budd- 
Vannimenus criterion jSj. 

3. Schrodinger equation and Hilbert space 

Difficulties of calculations of wave functions of the continuous spectrum are also present in 
problems of type of surface properties of metals where one has deal with the semi-infinite 
inhomogeneous electron gas. Attempts to replace the unbounded system by a system of finite 
size give rise to serious complications both in analytical and in numerical computations (see e.c. 
(§]). However, it is more instructive to analyze here the problem by an example of the tunneling 
in many-electron system. 

In this case the questions related to the eigenfunctions pertinent to the continuous spectrum 
of the Schrodinger equation ((J) can be considered with a fair degree of details and usefulness. In 
such a system the self-consistent Coulomb and exchange-correlation potentials effect inevitably 
on the shape of the tunnel barrier. The resultant contribution to tunnel current-voltage 
characteristics can vary from insignificant, as in metal-insulator-metal junctions, up to decisive, 
as in the Schottky-barrier metal-semiconductor junctions (see (Hj and references therein). The 
essential dependence of the transparency of self-consistent barrier on the energy of tunneling 
electrons and a reconstruction of the barrier with the applied bias voltage make in principal 
impossible the use of the tunnel Hamiltonian approach to describe such systems. Therefore it is 
necessary to formulate some regular scheme of tunnel current calculations that would form also 
a basis for the numerical realization of self-consistent solution. 

3.1. Scalar product and orthonormal basis in single-particle continuous spectrum 
Let two parts of system (left - L, right - R) occupy half-space z < (metal) and z > 
(semiconductor with a degenerate electron gas and the Schottky barrier), respectively. The 
effective potential energy V(z) of electrons is considered as independent of coordinates in the 
(x, y) interface plane. Then the single-particle Hamilton operator is: 



H = f + V(z), 



(20) 



where T is the kinetic energy operator. Let us assume following conditions for the asymptotes 
of V(z) 

V(z) -> V L , z -> -oo; V(z) -» 0, z -> oo. (21) 

In the case of metal-semiconductor junctions one can accept V(^) = V < for z < 0. Due 
to constraints 1)21 J) the eigenfunctions of continuous spectrum H have an oscillating behavior at 
the left or both infinities depending on the relation between the V L and the energy eigenvalue 
E. The second type eigenstates only give the contribution to the current. Let k > is the 
wavevector of wave function oscillations at z — > oo and q > is the same at z — > — oo. The 
eigenfunctions obey the equation 

HV E (x,y,z) =EV E (x,y,z) (22) 

and, in view of the translation invariance along the interface, they can be taken in the form 

ty E (x,y,z) = C fc Vfc(^)exp(ik||r||)/27r, (23) 

where E R = (k 2 + kn)/2 is the energy spectrum of electrons in the bulk of semiconductor. 

The wave functions of continuous spectrum should be normalized to the <5-function of quantum 
numbers. The Ea.(|23|) provides already the normalization to 5(k.\\ — kj,) in the lateral plane. 
This result (2ir in the denominator) is usually obtained by means of Born-Karman periodic 
boundary condition in the normalization box. Evidently, the system under consideration has no 
periodicity in z direction. To determine the constant C and to avoid cumbersome calculations 
of eigenfunctions for a finite-size system mentioned in Section^ let us define the scalar product 
of the eigenfunctions by 

/oo 
dzexp(-e\z\)ipl(z)ip kl (z) (24) 



with the natural definition of the eigenfunction norm by 

ll^fc 



2 = lim <V> fe h/> fel ) . (25) 

It is easily to check that the operator T in Eq. 1)201) and, therefore, Hamiltonian H are 
self-adjoined relative to the scalar product (|24jl. Hence, the orthogonality of eigenfunctions for 
different eigenvalues is provided by the self-adjointness of H. However, the quantum number 
k for current-carrying eigenstates is twofold degenerate. Thus any complex solution ijjf, and 
its conjugate form the linear-independent pair of solutions. In order to the degenerate pair 
of eigenfunctions do not violate the necessary orthonormality of the basis in the Hilbert space 
generated by Hamiltonian H they must be orthogonalized and normalized. 

Let us adopt as ip^ the wave function describing tunneling from right half-space to the left 
one and having asymptotes of form 

V£ = Cf (e- ikz + rfe ikz ) ,z - co; ^ = Cffi e"^, z - -oo. (26) 

The usual continuity condition of the probability flow density 

j k (z) = iP* k (z) (v + v+) Mz)/2 = const(z) (27) 

gives the relation between amplitudes of the transmission and reflection coefficients 



dE L (g,k ll ) 2 _ <9£ R (fc,k||) / 2 
dq " dk l 1_|rfc| 
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Here E L (q,k\\) and E R (k,k\ 
operator is denned by v = i 
momentum kn 



) are, respectively, the left and right energy spectrum, the velocity 
H, z . The conservation of the total energy E and the transverse 



E h (q,k {{ ) = E n (k,-k l{ ) = E (29) 
determines q(k) as a function k and vice versa. Accounting for Ea. (|29j) . the Ea. H28j) can be 



written in more compact form 



1 I**! 2 = i — kfc-l 2 - (30) 



dq/dk 

The contribution to the normalizing integral (|25j) is formed by infinite regions of the z axis 
where the asymptotic expressions (|26JI are valid. Using the definition (|24[) and the relation (|30j) 
the value of \Cf\ = 1/2tt can be found that provides 

Vm^f\^ 1 ) = d(k-k 1 ). (31) 

The second normalized solution ip^ which is linear-independent of can be obtained by the 
similar procedure with the following results 

Tpq = C\ (e^ + r^ z ) ,z^-oo; ^ = C^e ik \ z -> oo, (32) 
' ^| 2 = l-k a L | 2 (33) 



dk/dq 
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The other useful normalization 



C\ | = 1/27T, Urn «|< > = <5 (g - ft) , g = q(k) . (34) 



h(4)i4)) = i (*-*') ( 35 ) 

is obtained under condition 

, L|2 dg/dfc 

1^1 =^T~- ( 36 ) 

To elucidate a question about mutual orthogonality of ip^ and ^ needs to know the 

interrelation between the pairs (i]^, r^) and ^g(fe)> r g(fc)) ■ The required relations can be found in 
general form independently of a particular barrier if one realizes tyqnA as the linear combination 
of ipfc and that gives 

*g(fc) = dq/dk^' r 9( fc ) = _r ^* • ( 37 ) 



Using Eq. (|37|) one can show that terms like <5-function in scalar product {jPq^^ki/ can cel each 
other and hence 

The proof of Eq. (|38j) completes the construction of the orthonormal basis in Hilbert space of 
single-particle Hamiltonian H. 



3.2. Electron density and current density 

In the case of an equilibrium system the single-particle density matrix is p = p(H) and it 
has only diagonal non-zero elements p LL = F L and p RR = F R in the chosen basis due to the 
orthogonality relation 1)38)1 . Therefore, 

n(z) = Spph(z) = (39) 
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Since the Schottky barrier is completely situated in the semiconductor it is convenient to 
change the ^-integration in Eq. 1)40(1 by k- integration because k is the quantum number of the 
R-eigenstates. The result is 



n(z) = 2 dk 



dk,, {f l [^(fc,k||)] |^( 2 )| 2 + F R [JS? R (fc l k||)] |^)| 2 }. (41) 

Here the subindex q(k) in ip h was changed by k to recall the necessity to use the normalization 
(|55j)-(f37)|). 

The diagonal elements F L and F R of the density matrix determine respectively the occupation 
of L- and R-eigenstates and are described by Fermi distributions. The choice of the L- R-states 
as the basis allows to account mathematically for the presence of two independent reservoirs of 
particles (thermostats) in the left and right infinities. The Fermi level E F ' of each reservoir is 
determined by the proper neutrality condition far away from the interface and they are different 
as the bias voltage is applied to the junction. It is important to stress that both ifil{z) and 
ipu'(z) states extend to the both infinities and contribute to n(z) at any point. Thus the bias 
applied changes the position of each Fermi level. Substituting asymptotic expressions ()26)l and 
()32j) for tp R and in Eg. 1)41)1 and using the neutrality condition n(oo) = N R we obtain the 
equation determining the dependence of the Fermi level of electrons in semiconductor on the 
bias U = —V (compare with Eq.(3.101) in fTUj l 

A poo poo 

n ( z oo) = JV£ = ^ ^ dk J dk l{ F R [£ R (fc,k||)] + 
+ W?l™ dk /l dk l' ( 1_ ^ R ( fc ' k ll)] -^ R [^ R (^k||)]}- (42) 
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Here V is the structure voltage drop, k v = [2 max(£ L (0, 0) - £7,0)] , E F = Ef - U. At zero 
temperature E F = E R (kp) and this equation can be transformed into the equation for the Fermi 
wave vector of the right electrons 



kW* f dkdku / \ 

E(k,k\\)£[E F -U,E F ] 



The index R is suppressed here for short. It is easily to see from Eg. 1)43)1 that at U > we 
have kjy(U) > kp(0) and vice versa as it should be. Under the low barrier transparency (the 
reflection coefficient is of order of 1) the Eg .1)43)1 can be solved by simple iterative method 1 . 

1 In the case of metal-semiconductor junctions the relative correction to the solution of Eg. l43H due to a shift 
of the metal E]^ from the value — U by neutrality constraint is of order of » (U/E F )(1 — |rp;j 2 ) 2 and can be 
neglected in calculations of I — V characteristics of real structures with typical values 1 — |r^| 2 < 10 -4 , U < 1 eV, 
E F < 10 eV. Direct numeric calculations in Ref.|ll| have shown that a violation of the equality AE F — U 
becomes essential when the barrier height is < E F and barrier width is < 2n/k F . 



After the Fermi level is found one can calculate the tunnel current, averaging the current 
density (|27[) with density matrix and using the asymptotic representation of the wave functions 
at z — > oo. The result is 

'<">=-* r d n^ a i^-^^ ^ o-i-o- <«> 

Here F(E) is the Fermi distribution at T / 0, e = —1 is the electron charge and E(k,k\\) is the 
energy dispersion relation of the semiconductor electrons. 



4. Many-electron structures with finite number of electrons 

4-1- Many-electron atoms and molecules 

The neutrality problem of iteration procedure in self-consistent field theory has been considered 
here for the extended many-electron systems with undetermined number of electrons and the 
solution is suggested in Section |5J In the case of the finite many-electron systems like atom or 
molecule the number of electrons is exactly known and the insolubility of iteration equations 
seems not present. Anyway, the required charge state can be prepared if there is sufficient number 
of bounded-state levels. However, the question is then transformed into the impossibility to find 
the Poisson equation solution with a given angle symmetry if the charge density at the right- 
hand side does not have the desired symmetry. Fertig and Kohn |12j has considered the problem 
and discussed its source and consequences. But the solution suggested consists in additional 
artificial constraints on the sought variational solution. 

It seems, however, that separation of the induced charge in the self-consistent Poisson 
equation like it is done in Eq. ()13|) might also solve the problem in the case of finite Fermi 
systems. The induced charge can be taken in the form like in the Thomas-Fermi-Dirac theory 
[P5] with correlation potential included. 



4-2. Periodic electronic structures and supercells 

The spatial periodicity of crystal solids makes possible to replace numeric simulations of infinitely 
extended system by computations for a finite fragment (cell) supplemented by periodic boundary 
conditions. To map the preceding analysis of the interrelation between the charge distribution 
and boundary conditions for Coulomb potential on such a case let us rewrite the Eqs. ©,0 
for the finite interval (a, b) 

u (b) - u (a) = Air [ dzip(zi), (45) 



u(b) - it (a) = (6 - a)u' '(&) - An f dz\z\p{z\). (46) 

J a 

The periodicity demands the conditions for the total charge Q = J^dz\p(zi) = and the total 
dipole moment 

D = dz\Z\p{z\) = (b - a)u'(b)/An 



to be fulfilled. In the case of centrosymmetric structure there should be it' (a) = u'(b) = and 
therefore D = (see e.g. ^1] §§6, 13). At each ith iteration step it is easy to provide the 
neutrality condition = for the finite structure with a prescribed ionic charge by filling 
the necessary number of electronic states after solving the Schrodinger equation. However, the 
charge distribution with ^ is the rather likely result that cannot be prevented. In this 
case the periodic boundary condition for the solution of the Poisson equation can be fulfilled 
only at u' ^ at the boundaries. At the same time the behavior of the potential should be 



essential contorted in comparison with a "true" charge distribution at D = that produces 
increasingly distorted charge distribution at the next iteration cycle. In the case of Fourier- 
series decomposition of the electron density and potential, the discontinuity of the potential at 
the boundaries entails a growth of the solution u and, therefore, its derivative v! again near the 
boundaries owing to the Gibbs effect |15| . 

In the case of two- or three-dimensional cell the Eqs. (j45j) and (|46|) should be replaced by 

j) dSv s -Vu = 4tt J drp(r) = 47rQ cc u, the Gauss theorem, (47) 

"Sccll ^ccll 

j> dS[v s -Vu-u(v s -V)]Y = 4vr J drrp(r) = 47rD ccU , (48) 

^>ccll ^ccll 

with the same conclusions as in the one-dimensional case above. Here O ce n and S ce \\ are the 
volume and bounding surface of the cell, respectively, v s is the external unit normal to the cell 
surface. The formulae (|47[) and (|48|) are derived from the Poisson equation with the use of the 
Green theorem 

Au(r) = 47rp(r), / dv^Ml = / dr^AQ+J) dS ($W - W$) (49) 
Jn Jn Jsn 

that requires for proper differentiability of the involved functions and correspondent smoothness 
of the bounding surface S to be fulfilled [T5] . 

The incompatibility of the boundary conditions with the right-hand-side of the Poisson 
equation makes the Dirichlet problem as ill-posed one with corresponding strong perturbance 
of the solution in response to insignificant perturbance of the initial data ^7]. This mechanism 
may be responsible for the observed deterioration of current self-consistent iterative algorithms 
manifested in the long- wavelength charge instability, which is named as "charge sloshing" |18j - 

It is necessary to note that at initial iteration step the solution of self-consistent Poisson 
equation (|T3|) with the induced electron distribution n- m( ± defined by Eq. (|TT|) should give a good 
starting guess for the potential and the valence electron distribution in the cell since they are 
self-consistent and meet the boundary conditions. However, the expression (|11|) is evidently 
inappropriate in the case of crystals with the filled energy bands and it needs to use a next quasi- 
classical approximation for the induced electron density expressed as a function of derivatives 
of the potential (see, e.g. [2Tj ) . 

Let us consider now the application of the results of Sec. El to the continuum spectrum of 
infinitely extended periodic system. The eigenfunctions of single-particle Hamiltonian according 
to Bloch theorem can be taken in the form 

^j k (r) = Cexp(ikr)</> jk (r), (50) 

where ^jk( r ) is the cell-periodic part, j is the energy band index, k is the wave vector, and C is 
the normalizing constant. The definition of the scalar product like in Ea. (|24j) by the expression 

(^jklV'jikx) = Hm j drexp(-e|r|)^ k (r)^ jlkl (r) (51) 

and the eigenfunction norm like in Eq.((2SJ) provides the self-adjointness of single-particle 
Hamiltonian and results in the orthonormal basis of the Hilbert space with \C\ 2 = (2n)~ 3 



(V'jkl^hki) = <5 (k — k x ) <5j jl 



(52) 



under natural conditions k, ki € 1st Brillouin zone and 



J dr |^ jk (r)| 2 = Q cell . (53) 

f^cell 

The relationships just obtained allow to determine the electron density inside the cell by 

n W = 2 E / ^hk(r)| 2 e(£ F -£j(k)), (54) 

where the multiplier 2 takes into account the spin degeneration and ^-integration is taken over 
the unit cell of the reciprocal lattice that is the first Brillouin zone. Because the volume of the 
Brillouin zone 13 BZ = (27r) 3 /O ce u we have 

f dk 1 

J (2^)3 "IW (55) 

Thus the band filling and the position of the Fermi level can be determined from the neutrality 
condition 

J drn(r) = N+, (56) 

^cell 

where N + is the ion charge of the unit cell. 

The expressions (|53|) -(|56 |) allow to abandon the use of unnecessary Born-Karman boundary 
condition for the eigenfunctions V ; jk( r )> which limits admissible points in /c-space by the discrete 
set. As a result, one can calculate the energy bands £j(k) of perfect crystal or make integrations 
over Brillouin zone, using any set of /c-points dictated by selected algorithm [22], and avoid 
disadvantages of the cell size extension beyond the size of the minimal unit cell to increase the 
sampling density of /c-points [23] . 

Similarly, the artificial periodicity of the supercell JH] can be eliminated from the calculations 
of solid surface or imperfect crystal with a point defect and replaced by asymptotic boundary 
conditions at infinity for the Coulomb potential and wave functions. This removes the known 
difficulties introduced in calculations of extended systems by making use the slab [21] or supercell 
(2H]-!2S] geometries. In this case the induced charge resulted from the free carriers or nonuniform 
polarization of the valence electrons should be introduced into the self-consistent Poisson 
equation using the effective mass approximation or the macroscopic electric susceptibility, 
respectively. 

5. Concluding remarks 

The origin of all difficulties with self-consistency listed above is the long-range character of the 
Coulomb interaction that is responsible for an interdependence of the charge distribution with 
the boundary conditions on the potential. Thus large-scale self-consistent distributions can only 
result from the direct solution of the self-consistent Poisson equation itself because the boundary 
conditions take into account the reaction of distant charges that are exterior ones relative to 
the considered system. The separation of the induced charge as a function of the potential 
and the modification of the Poisson equation is only way in order to obtain the self-consistent 
distributions of potential and charge with due accounting for the boundary conditions. 

To elucidate the point in more detail let us consider the electrostatic part of the Kohn-Sham 
energy functional 

*.-\f ^r'f^l (57) 



and the corresponding contribution to the effective potential in the Schrodinger equation |22j 



dr'J^-. (58) 
|r — r| 

At first glance, there is no need for the Poisson equation since Eqs. (|57[) and (|58[) give us already 
the explicit relationships between the necessary quantities and the charge density. However, the 
function 1/ |r — r'| in the integrand of the Eq. (|58f) is the Green function of the Laplace equation 
with the zero boundary condition at the infinity. The expression (|58|) is the true solution of the 
Poisson equation if there are no charges outside the integration region. Evidently, this is not 
the case for infinitely extended systems or when periodic boundary conditions are specified. 
Let us assume that there are two subsystems with non-intersected charge distributions 

p(r) = Pl (r) + p 2 (r); Pl (r) + 0, r G T x ; p 2 (r) ± 0, r G Y 2 ; Yx n Y 2 = (59) 

and, respectively, 

cp(r) = cpi(r) + tp 2 (r); Aip lj2 = -4vrpi i2 , (60) 

where ip is the Coulomb potential. Substituting the expression (|59|) for p in Eq. 1)5 7 Jl and using 
Poisson equations (|60|) together with Green theorem (|49|) we obtain 

E es = J dvpiLp + J dripAip + / drp 2 <p + J dnpAip = E cs i + E es2 . (61) 

Ti Ti T2 T2 

Now it is possible to consider the total energy -Etoti of the subsystem 1 as the functional of 
/Oi(r), ^fEi( r ), and <p(r). Then the necessary conditions of the functional minimum are [2Z1 

5£W**li ( r ) = 0, ^esi/^(r) = (62) 

that must be supplemented by the expression for the effective potential u e fji = 5E es i/5pi + n xc i. 
The first equality gives rise to the Schrodinger equation and the second one is just the Poisson 
equation. The variation Sip and the sought potential tp have to satisfy such boundary conditions 
that ensure 

j> dS (<pV5ip - SipVcp) = (63) 

at the surface confining the region Ti only. It implies that boundary conditions to Poisson 
equation must be properly fixed during the iterative process. If a realization of the equality 1)63(1 
at given boundary conditions turns out to be impossible then the separate investigation of two 
subsystems is incorrect. This derivation shows that the solution of the Poisson equation cannot 
be replaced by direct variations of Eq. (|57|) - (|58J) in the course of searching for the E tot i minimum 
for the infinitely extended electronic systems with the periodic or assigned asymptotically at 
infinity boundary conditions. 

It is desirable to add some comment on the widely- used "mixing" method of fight against the 
charge instability. The very essence of the method lies in the use of some linear combination of 
results of previous iterative steps to make up the input for the next step. Such an approach is 
named as "iteration with memory" in the iterative calculus and it is destined to accelerate 
the convergence rate of the iterative process |28| . But this procedure cannot transform a 
divergent iteration scheme to convergent one. There are many reasons, including mentioned 
in the present work, to believe that charge instability observed in simulations of large electronic 
structures is rather sign of divergency than slow convergency of the simple iteration cycle 
of Poisson^Schrbdinger— >Poisson steps. The continuous appearance of new more and more 



cumbersome and sophisticated mixing methods during the last three decades is the best evidence 
of this point of view (see some historic commentary in pQ, |18|-|2f)j. |29j). As a rule, the 
demonstration of the benefit of a new method is accompanied by illustrations of the lack of 
convergence of the old one as the structure size becomes larger. It is of interest, the mixing 
scheme based on some handmade modeled screening succeeds relatively [22], although it was 
constructed on the ground of rather formal mathematical reasoning than the underlying physics 
discussed here. 

It needs also to note that mixing schemes may produce a spurious convergency (see |2Uj . 
p. 11176). Thus it is necessary to check whether the norm of the functional derivatives 
||5£ , cs i/5(/3(r)|| and ||<5i?toti/^j5;i(r)||, which are residuals of Poisson and Schrodinger equations, 
are minimal along with the total energy i^toti- Of course, i?toti is the total energy of sufficiently 
large but finite system that can be approximated by assignment of boundary conditions as in 
the infinitely extended one. 
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